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Abstract 
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1 Introduction and motivation 

1.1 The abstract formulation of the problem 

Throughout this paper, we assume that n G {2,3, . . .} and that 




We also assume we are given n strictly increasing breakpoints on the real line: 




Our goal is to 

(3a) find a vector x = (xi, . . . , x„) G X 

such that 

(3b) (ii,xi),...,(i n, Xn) satisfies a given set of constraints. 

Note that for every x G X, the pair {t,x) G X x X induces a corresponding piecewise 
linear function or linear spline (see ||27|| and ||47|) 

T — t- 

(4) fn ^) : [h, tn] ^ 'R: T ^ Xi + {Xi+i - Xi)- 

where t G [i;, ij+i] and i G {1, ...,n — 1}, which passes through the points 

{{U,x^) elR^\iel}. 



The set of constraints mentioned in (|3b|) will involve the function f(^t,x)- Let us list sev- 



eral types of constraints which are motivated in Section [L2l below: 

• interpolation constraints: For a given subset I of {1, . . . ,n}, the entries X/ are pre- 
scribed: (Vz G I) f(^tx)iU) = Vi- This is an interpolation problem for the points 
{(i„y,) G1R2 I 

• slope constraints: For a given subset J of {1, . . . , n — 1} and for every i G I, the sZope 
(5) s,- := ^i±l^ 

of /(f I [t;,t;^j] must lie in a given subset of R. 



• curvature constraints: For a given subset I of {1, . . . ,n — 2} and for every i G I, 
|s/_|_i — Si\, the distance between the slopes s, and s,+i of two adjacent intervals 
[tj, tjj^i] and [tj+i, ^;+2]/ rnust lie in a given subset of R. 

1.2 A concrete instance in road design 

The problem © introduced above and its solutions have several direct applications in 
engineering and computer-assisted design. For instance, an engineer may want to verify 
the feasibility of a design, or adapt the design according to the constraints. Examples 
drawn from Computer- Assisted Design (CAD) include designs for roadway profiles, pipe 
layouts, fuel lines in automotive designs such as cars and airplanes, overhead power 
lines, chairlifts, cable cars, and duct networks. 

Our primary motivation for this work is automated design of road alignments. A road 
alignment is represented by the centerline of the road, which is idealized as a (generally) 
nonlinear, smooth curve in R"^. To facilitate construction drawings, civil engineers reduce 
the three-dimensional road design to two two-dimensional parts, horizontal and vertical. 

The horizontal alignment is the plan (or map) view of the road. In the vertical view, the 
ground profile g: [ti, f„] — > R shows the elevation values of the existing ground along the 
centerline (see the brown curve in Figure [T]). Since earthwork operations such as cut and 
fills are expensive items in road construction, a civil engineer tries to find a road profile 
represented by a linear spline f(^t,x)) that follows g as closely as possible. 
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Figure 1: Vertical profiles of ground profile g, initial road design f{^t,v)r ^rid final road 
design f{^t,w) ^r a highway design near Kelowna, B.C., Canada. The solution was found 
with the ParDyk algorithm (Algorithm |4.5l below) using a design speed of 130 km/h, and 
a maximum slope of 4%. (These engineering constraints translate into a specific set of 
slope and interpolation constraints.) 
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Design constraints imposed on /(t^x) by the engineer or by civil design standards such 
as those published by the American Association of State Highway and Transportation Officials 
(AASHTO) 111 include the following: 

• At a certain station ti, the engineer may have to fix the elevation to allow for 
construction of an intersection with an existing road that crosses the new road at t^. 
This corresponds to the (mathematically affine) interpolation constraint mentioned 
in Section im 

• For safety reasons and to ensure good traffic flow, AASHTO requires that the slope 
between two stations ti and is bounded above and below. These are the (math- 
ematically convex) slope constraints of Section [Lll 

• In a road profile, engineers fit a (usually parabolic) curve at the point of intersection 
of two line segments. The curvature depends on the grade change of the line seg- 
ments and influences the vertical acceleration of a vehicle, as well as the stopping 
sight distance. AASHTO requires bounds on the slope change. This corresponds to 
the (mathematically convex) curvature constraint in Section [Lll 

• In some cases, the engineer requires a minimum drainage grade to allow flow and 
avoid catchment of storm water. These (mathematically challenging nonconvex) 
slope constraints are discussed in Section |5T] below. 

We denote the starting spline for a road profile (see the cyan curve in Figure [T]) by 

(6) where f = (fi,...,f„) ^ ^ and u = (z7i,...,z;„) G X. 

In practice, the starting spline could simply be the connected line segments for the inter- 
polation constraint, or it could be generated from the ground profile g by using Bellman's 
method (see ||11| and ||45| for details). In either case, we assume that t is given and fixed, 
and that we need to decide whether or not /(^ is feasible with respect to the aforemen- 
tioned constraints. If v leads to a feasible spline, then we are done; otherwise, we wish to 
find zf G X such that the new road spline jy) (see the blue curve in Figured) satisfies 
the design constraints. Ideally, f{^t,w) is close to the ground profile represented by {t,v). 
Finally, if there is no w G X making the problem feasible, then we would like to detect 
this through some suitable measure. 

1.3 Main results and organization of the paper 

We now summarize the main contributions of this paper. 
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In principle, there are numerous constraints to deal with for the problem ^ in the 
context of road design. Fortunately, the constraints have a lot of structure we can 
take advantage of and we demonstrate that the constraints parallelize which allows to 
reduce the problem to six constraint sets (see Section IZS]) each of which admits a closed 
form projection formula (see (O, ([21]), and (|31]), (|99l) ). 



We provide a selection of state-of-the-art projection methods, superiorization algorithms, 
and best approximation algorithms (see Sections |3] and S]), and adapt them to the road 
design problem. 

We present various observations on the algorithms and their relationships (see Re- 
mark |331 Remark 1121 Remark I4l5l Remark l4T6l and Example O) 

We report on broad numerical experiments (see Section (6]) introducing for the first 
time performance profiles for projection methods. 

Based on the numerical experiments, we recommend CycP+, which is an intrepid 
form of the method of cyclic projections, as an overall good choice for solving feasi- 
bility and best approximation problems. 



The remainder of the paper is organized as follows. Section |2] contains a detailed anal- 
ysis of the projection operators encountered in the road design problem. We take ad- 
vantage of aggregating constraints and derive simple formulas. Projection methods for 
feasibility problems are reviewed in Section |3l Because we are working with more than 
two constraint sets, we adapt to this situation by working in a product space if needed. 
If more than just a feasible road design is desired, then the engineer has to consider op- 
timization algorithms. We review two types of such methods (superiorization and best 
approximation algorithms) and adapt them to our problem in Section 13.51 and Section H] 
respectively. Nonconvex constraints are investigated in Section |5l We report on numerical 
experiments in Section [6l The final Section [7| concludes the paper. 

For notation and general references on the mathematics underlying projection methods, 
we refer the reader to the books M, El/ BH/ IH, IH, and BH. 



2 Constraints and projection operators 
2.1 The projection onto a general convex set 

In this section, we make the constraints encountered in road design mathematically pre- 
cise. Almost all of these constraints turn out to lead to sets that are convex and closed. 
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Recall that a set C is convex if it contains all line segments between each pair taken from 
C: 

(7) (Vco G C) (Vci G C) (VA G [0, 1] ) (1 - A)co + Aci G C. 

If C is a nonempty closed convex subset of X, then for every x G X, the optimization 
problem 

(8) C) := min ||x — c||, 

ceC 

which concerns the computation of the distance d{x,C) from x to the set C, has a unique 
solution, denoted by Pqx and called the projection of x onto C. The vector P^^ is charac- 
terized by two properties, namely 

(9) PcX G C and (Vc G C) (c - PcX,x- Pqx) < 0. 

(For a proof, see, e.g., HI Theorem 3.14]). The induced operator Pc ^ ^ ~^ C is the projec- 
tion operator or projector onto C. There are several examples that allow us to write down 
the projector in closed form or to approximate Pq provided C is the intersection of finitely 
many simple closed convex sets (see, e.g., HJ Chapters 28 and 29].) In the road design 
application, it is fortunately possible to obtain explicit formulas; in the following subsec- 
tions, we will make these formulas as convenient as possible for software implementa- 
tion. Let us do this for each of the three types of constraints. 



2.2 Interpolation constraints 

We assume that J is a set such that 

(10) {l,n} C 7 C {l,2,...,n}, and y = (y,) g 
is given. Set 

(11) y := |x = (xi, ...,Xn) G X I (Vz G I) Xj = yi}. 

The closed set Y is an affine subspace, i.e., (Vy G Y) (Vz G Y) (VA G R) (1 — A)y + Az G Y; in 
particular, Y is convex. For convenience, we record the explicit formula for the projection 
onto Y. 

Proposition 2.1 (interpolation constraint projector) The projector onto Y is given by 
(12a) Py:X^X 

'y/, ifi G I; 



(12b) {Xi, X2, . . . , Xn) ^ (Ci , C2, ...,Cn), where c, 



Xi, ifi G {1,2,. . .,n} \ J. 



6 



2,3 Slope constraints 



2.3.1 A useful special case 



We start with a simple special case that will also be useful in handling our general slope 
constraints. To this end, let z G {1,2, . . . , n — 1}. The constraint set S, imposes that the 
absolute value of the slope /(f^) for the interval [ti, tj^i] is bounded above, i.e., there exists 
a;, > such that 



(13) 



Si := [x = (xi, . . . ,x„) G R" I \xi+i — Xi\ < OLi}. 



Indeed, if the actual maximum absolute slope is cr, > 0, then, setting a/ = cri\tij^i — ti\, we 
see that 



(14a) 

(14b) 
(14c) 



U+1 — U 



< CTi 4^ |Xj+i — Xi\ < Oii ^ ±(^f+l — ^i) < Oi-i- 

^ ± (Cf+i - Ci, X) < Xi 

<^ —Oii < — Ci, x) < Oii, 



where and denote the standard unit vectors in X (with the number 1 in position i 
and i + 1, respectively, and zeros elsewhere). The last characterization reveals that is the 
intersection of two halfspaces whose boundary hyperplanes are orthogonal to — e/. 
In particular, S, is a closed convex subset of X. Using e.g., [6, Example 28.17], we obtain 
for every x G X, 



(15) 



X + 
x + 



-Oii 



{ei+i - Ci, X 



- ei), if - e,,x) < -a;; 



if —Oii ^ ~ ^iiX) ^ 



This formula shows that (Ps,^)/ = ^7 for every; G J \ {i, z + 1}. Thus, the only entries that 
possibly change after executing P5. are in positions i and z ' + 1; after some simplification, 
we obtain for these entries the formula 

{\{xi + X/+1 + Oii, Xi + ^z+i - Oii), if Xi - X/+1 > Oii; 
(x;,x/+i), if |x/+i -Xi\< Oii) 

\{Xi + Xf + l - ft;, Xi + X; + l + a,), if X; + l " X, > Oii- 

We note that 

(17) x^S, ^ (Ps,x), 7^ X, and (Ps 

furthermore, if a, = +00, i.e., no slope constraint, then (|T6)) is valid as well. 



2.3.2 The general case 



Now we turn to the general case. We assume the existence of a vector a = [di) G R" 
such that the constraint set is 

(18) S= fl S, = {(xi,...,x„) gR" I (VzG {l,...,n-l}) |x,+i-x,| <a,}, 

ie{l n-l} 

where Si is defined in (|T3|) . While we obtained an explicit formula to deal with a single 
slope constraint (see ((T6|) ) we are unaware of a corresponding formula for P5. Further- 
more, since Pg- possibly modifies the vector in positions i and i + 1 (but not elsewhere), 
we cannot use (|T6|) for the sets Sf and S;_|_i concurrently because their projections possibly 
modify positions {i, i + 1) and {i + 1, z + 2) (see ^7\ ), but not necessarily in a consistent 
manner at position z ' + 1! However, by combining the n — 1 slope constraints according 
to parity of indices, i.e., by setting 

(19) Seven := Q S;, and Sodd := H 

i6{l,...,n-l}n(2N) !€{!,.. .,w-l}n(l+2N) 

we see that 



(20) S — Seven H Sodd 

can be written of the intersection of just two constraint sets! Furthermore, (fT6]) yields the 
fully parallel update formulas: 

Proposition 2.2 (convex slope constraint projector) For every x G X, the projectors onto 
and Sodd i^^s given by 

(21a) PSeven^ = i^l' (^^52^)2, (^52^)3, (^^54:^)4, (^54:^)5, • • • ) ^ X, 



where the last entry in (|21a|) is %„ ifn is even; and 

(21b) PS,,,X= ((PsiX)i,(PsiX)2,(Ps3x)3,(Ps3x)4,...) G X, 

where the last entry in (|21b|) z's x„ z/n zs odrf. 



The constraints making up the aggregated slope constraint are very special polyhedra, 
namely "stripes", i.e., the intersection of two half spaces with opposing normal vectors. 
This motivates the technique, which originates in [42J (see also [16L [ilZJ, and [41J), of not 
just projecting onto these sets but rather inside them: either we reflect into the strip or (if 
we are too distant from the strip) we jump to the corresponding midpoint of the strip. Let 
us record the formula for this operator. 
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Proposition 2,3 (intrepid slope constraint projectors) The intrepid counterpart of (|T6)), i^ 

(22) ixi,Xi+i) <j {xi,Xi+i), if\xi+i -Xi\ < oca 

{xi+i,Xi) + sgn(x/+i - Xi){-ai,ai), if cti < \xi- < Idi. 



These operators lead to the intrepid counterpart of (|2T|). 

2.4 Curvature constraints 
2.4.1 A useful special case 

Again, let us start with a simple special case that will also be useful in handling our 
general curvature constraints. To this end, let z G {1, . . . , n — 2}. The constraint requires 
that the difference of consecutive slopes is bounded above and below. Hence there exists 
7, G R and 5i G R such that 7/ > 5i and 

(23) Q:=|x=(xi,...,x„)GR« 7,. > ^i±l^^ - ^i±i^ > 
Set 

(24) Ti := if+i - ti, := i/+2 - ^z+i/ and := T^+ie, - (t, + T,+i)e,+i + T;e,+2- 
Then > 0, T;+i > 0, and for every x G X, 

(25) X G Q 7iTf+iT, > (uf,x) > SiTi+iXi. 

Again, we see that C; is a closed convex polyhedron and by e.g. ||6l Example 28.17], we 
obtain 



(26) Pc,::^^<^ 



X, if (5,T,T;+i < (w,-,x) < 7zTiT;+i; 

7,T,T; . 1 - (Uf, X) . . , , 



Since G span{e;,e;4_i, ef_|_2}, it follows that 

(27) {;■ G {l,...,n} I x^- ^ (^'c,^);} ^ {2,2 + 1,2 + 2}. 



Here sgn denotes the signum function. 
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2.4.2 The general case 



Now we turn to the general curvature constraints. We assume the existence of a vector 
c = {ji) G R"~^ and d = (Si) G R"~^ such that the constraint set is 

(28a) C = fl Q 

ie{l,...,n-2} 



(28b) 



I ti+2 — — ti J 



Because of (|27l) , we can and do aggregate these n — 2 constraints into three constraint sets 
that allow projections in closed form. To this end, we set 

(29) (V/G {1,2,3}) C[^.]:= f) Q 

!6{l,...,«-2}n(/+3N) 

so that 

(30) c = C[i]nq2]nq3]. 



Combined with (|26|) , we obtain the following. 



Proposition 2.4 (curvature constraint projector) For every x G X, i/ze projectors onto C^^j, 
C[2], and C[3] are given 

(31a) Pc[j,X = ((Pci^)i, {Pc,x)2, (Pci^)3, (PC4^)4, (PC4^)5, (PC4^)6, • • . )' 

(31b) PcpjX = (Xi, {Pc,x)2, (PC2X)3, (PC2X)4, (Pc5^)5, (^€5^)6, (^€5^)7, • • • )' 

(31c) PcpjX = (xi, X2, (PC3^)3/ (^^€3^)4, (J'C3^)5/ {PCeX)6, (^'Cg^)?, (f'C6^)8/ • • • )' 

respectively, where each Pq. is given by (|26l) . 



We now record the intrepid counterpart. 
Proposition 2.5 (intrepid curvature constraint projectors) The intrepid counterpart of (|26l) 
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IS 



X - „ — — Ui, if {u„ x) < TiTi+i {33i - ji) /2; 



(32) X ^ < 



,.\\2 



^ {Uj, X) - ^/T,T,+i 



ifTiTi+i{3Si-ji)/2 < {ui,x) < SiTiTi+i; 
ifjiriTi+i < {Ui,x) < T,-T,-+i(37/-<5,-)/2; 



These counter-parts induce the intrepid variant of (|3T]). 



2.5 The convex feasibility problem 

The convex feasibility problem motivated by road design is to find a point in 



(33) 



C :— y n Seven H Sodd ^ C[2] H C^j, 



where the sets on the right side are defined in (|TT|) , ((T9|) , and (|29l) . Note that we have 
explicit formulas available for all six projectors (see Propositions 12 . 1 112 .21 and I2.4[) . More- 
over, we have intrepid variants of the projectors onto the fiv^ constraint sets Seven/ Sgdd/ 
C[x], Cp], and Cpj (see Propositions |2.3l and|2.5 [) . 



3 Feasibility problems and projection methods 



In this section, we present a selection of classic projection methods for solving feasibility 
problems. To this end, let m G {2,3, . . .} and C\, . . . ,Cm be nonempty closed convex 
subsets of X. 



The intrepid variant of Py is Py itself because the interior of the affine subspace Y is empty. 



11 



3.1 The m-set feasibility problem and its reduction to two sets 



Our aim is to 

(34) find x e C := CiH ■ ■ ■ nCm 7^ 0, 
or, even more ambitiously, to compute 

(35) Pcv 

for some given point u G X, or some intermediate between feasibility and best approx- 
imation. Of course, we have the concrete scenario (|33l) of Section |2.5[ where m = 6, in 
mind; however, the discussion in this section is not limited to that particular instance. 



Projection methods solve (|34|) by generating a sequence of vectors by using the projec- 
tors Pcj, . . . , Pc„, onto the individual sets. (For further background material on projection 
methods, see, e.g., dl, |[T5|, and ||22|.) As is illustrated by Section |2l these projectors are 
available for a variety of constraints appearing in practical applications. Before we pro- 
ceed to catalogue the algorithms, we note that some algorithms (see, e.g.. Algorithm 13.121 
below) work intrinsically only with two constraint sets. Because m in our application is 
not large, this turns out to be not a handicap at all as we can reformulate (|34l) in a product 
space in the following fashion: In the Hilbert product space 

(36) X := X'«, 

equipped with 



(37) (x,y) = ^(x„y,) and ||x|| = "^11-112 

i=l 

where x = (xi, . . . ,Xm) and y = (1/1, . . . ,ym) belong to X, we consider the Cartesian product 
of the constraints together with the diagonal in X, i.e., 

(38) C:=CiX---xCm and D := {(x,. . G X | x G X}. 
Then for every x G X, we have the key equivalence 

(39) xgC ^ (x,...,x) G CnD; 

thus, the original m-set feasibility problem (|34l) in X reduces to two-set feasibility problem 
for the sets C and D in X. This approach is viable because the projectors onto C and D 
are given explicitly (see, e.g., [6^ Propositions 28.3 and 28.13]) by 

(40a) Pq ■■ (xi, . . . , Xyn) ^ {Pci^i, PCn^m) 

and by 

(40b) Px)-. {xi,...,Xm)^ {y,...,y), where y = ^(xi H hx,«), 

respectively. 
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3.2 Projection methods and . . . swiss army knives! 

Projection methods use the projectors onto the given constraint sets in some fashion. Be- 
cause the squared distance functions, d^(-,Q), are Frechet differentiable with derivative 
2(Id— Pq) (see, e.g., [|6> Corollary 12.30]) projection methods are first-order methods — 
these are known to allow at best for linear convergence. Not surprisingly, they are not 
always competitive with special-purpose solvers ||37j| ; however, when projection meth- 
ods succeed (see 1161 for a compelling set of examples), then they have a variety of very 
attractive features: 

• Projection methods are easy to understand. This is important an industry, where 
mathematical/algorithmic considerations are only one part of an engineer's job. 
The engineer will not typically be familiar with the latest research developments in 
all branches of relevant mathematics. On the other hand, the idea of a projection 
method is often immediately grasped by drawing some sketches. 

• Projection methods are easy to implement. In Section |2l we have seen various formu- 
las for projection methods. These involve simple operations from linear algebra and 
are easily programmed. 

• Projection methods are easy to maintain. Once implemented, the code for these meth- 
ods is typically small and in fact smaller than other pieces of code dealing with data 
input/ output. This makes maintenance straightforward. 

• Projection methods are easy to deploy. Because of typically small memory require- 
ments, it makes them much easier to deploy on low memory computers like mo- 
bile devices. Also, the code base required for projection methods is in most cases 
significantly smaller than the size of libraries for linear or nonlinear optimization 
solver software. Thus, projection methods satisfy some of the key requirements for 
embedded optimization [13J, where the solution of one method is used within an en- 
compassing algorithm. 

• Projection methods are inexpensive to implement. Because of typically straight for- 
ward implementations, there is no need for commercial optimization solver soft- 
ware. 

• Projection methods can be very fast. If the iterations in a projection method can 
be executed quickly, then for certain classes of problems projection methods can 
become very competitive with traditional optimization algorithms. In Section 16.61 
below, we illustrate the enormous potential of projection methods when compared 
to algorithms for linear programming or even mixed-integer linear programming. 
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In summary, projection methods possess the same essential characteristics ofswiss army knives: 
they are flexible, lightweight, simple and very convenient provided they get the job done. 
If the saw included with the swiss army knife in your pocket cuts the branch of the tree, 
then there is no need for you to either get a big saw out of the garage or to buy a chainsaw 
from the hardware store! The road design feasibility problem analyzed in this paper adds 
a new compelling success story. 

3.3 A catalogue of projection methods 

In this subsection, we provide a list of projection methods. Each of these methods pro- 
duces a sequence that converges — sometimes after applying a suitable operator — to a 
point in C provided that C 7^ (and perhaps an additional assumption is satisfied). Basic 
references are the books lU, ||T5l, 1221, [281, 13511 , and 13611 ; if these books are insufficient, 
we include further pointers in accompanying remarks. 

Note that numerous generalizations of projection methods are known. These typically 
involve additional parameters (e.g., weights and relaxation parameters). Because we shall 
compare these methods numerically, we have to restrict our attention to the most basic 
instance of each method. 

While these methods are generally not able to find Pcv, they may lead to feasible solu- 
tions that are "fairly close" to v provided the starting point is chosen to be v. 

We start with the method of cyclic projections, which has a long history (see, e.g., 1128 II ). 

Algorithm 3.1 (cyclic projections (CycP)) Set xq = v and update 

(41) (VfcGN) Xfc+i:=rxfc, where T := Pc,„Pc,„_, ■ ■ ■ Pc,- 

A modem variant of CycP replaces the projectors in (|4T1) by intrepid counterparts when 
available. 

Algorithm 3.2 (cyclic intrepid projections (CycP+)) Set xq = v and update 

(42) (Vfc G N) Xfc+i := Txk, where T := RmRm-i ■ ■ ■ Ri, 
where each Ri is Pq or an intrepid counterpart thereof. 

Remark 3.3 (CycP+) CycP+ (also known as ART3+) is available for (l33l) because we have 
the intrepid counterparts of the projectors at our disposal (see Propositions 12.31 and 12.51 
and Footnote l2|). This method is fairly recent; see [|16B , (121, and [,42 , 1 (for constructions 
based on halfspaces and subgradient projectors), and also H (where one of the sets is an 
obtuse cone K so that the intrepid projector is actually the reflector 2Pj: — Id). 
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While Algorithm 13. II is sequential, the following method is parallel: 
Algorithm 3.4 (parallel projections (ParP)) Set xq = v and update 

(43) (VfcGN) x^+,:=Tx,, where T := l(Pci + • • • + PcJ- 

Remark 3.5 (ParP) In view of (|40|) , it is interesting to note that ParP is equivalent to iter- 
ating Pd^*C/ i-G-, to applying CycP to the subsets C and D of X. See also [46, Corollary 2.6] 
and 121 Section 6]. 

The next method can be seen as a hybrid of CycP and ParP. 
Algorithm 3.6 (string-averaging projections (SaP)) Set Xq = v and update 

(44) (VfcGN) Xk+i:=Txk, where T := 1 (P^ + + • • • + Pc„, • • • i'cAi)- 
Remark 3.7 (SaP) For further information, see | [20| and [|2T| (and also ||3l. Example 2.14]). 

Algorithm 3.8 (extrapolated parallel projections (ExParP)) Set xq = v and update, for 

every G N, 
(45) 



Xj^j^i := Txj^, where Tx :- 




sr^m 1 

2^i=l \ 


\x — Pci^W 


2 




[x - Pc,x) 


2 



^(PqX — x), if X ^ C; 

i=l 

otherwise. 



Remark 3.9 (ExParP) This method is actually an instance of the subgradient projection 
method applied to the function x i— ^ J2lLid^{x>Ci); see |24| for further information. 

Algorithm 3.10 (extrapolated alternating projections (ExAltP)) Assume that Ci is an 
affine subspace — this is the case in our application when we choose the interpolation 
constraint. Set xq = v, and let G N. Let 4 be a nonempty subset of {1, . . . , m} contain- 
ing exactly indices such that for each ; G {2, . . . ,m}, j belongs to Jjt frequently. Given 
Xfc, set 
(46a) 

^^^^ I 1, otherwise, 

and then update 

(46b) Xfc+i = z/c + jikiVk - Zk)- 
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Remark 3.11 (ExAltP) See (Zl Algorithm 3.5] for further information on this method. In 
our implementation, we chose Ijt = {2, . . . ,m} so that mj^ = m — 1. 



The next method is different from the previous ones in two aspects: First, it truly op- 
erates in X and thus has increased storage requirements — fortunately, m is small in our 
application so this is of no concern when dealing with (|33|) . Second, the sequence of in- 
terest is actually different and derived from another sequence that governs the iteration. 

Algorithm 3.12 (Douglas-Rachford (D-R)) Set xq = (v, . . . ,v) e X = X"" . Given G N 
and Xfc = . . . , e X, update to x^+i = (x/c+i,i, . . . , x^+i^^), where 



^ m 

(47a) Xk= — V] Xk^i 



m 

i=l 



and 

(47b) (Vz G {1, . . . , m}) Xfc+i^f = x^,; -Xk + (2xfc - xj,j). 

The sequence of interest is not (xjt)fc6N but rather (xfc)fce]N. 

Remark 3.13 (general Douglas-Rachford algorithm) Let us briefly sketch how the up- 
date formula (|47|) is derived from the general splitting version of D-R, which aims to min- 
imize the sum / + g of proper lower semicontinuous convex functions / : X — > ] — oo, +oo] 
and g: X — ] — oo, +oo] . Given Xq G X, the algorithm proceeds via 

(48) (Vfc G N) Xk+i = Txk, where T = Prox^(2 Prox/ - Id) + Id - Prox^, 

and where Proxy denotes the proximal mapping (or proximity operator) of /, i.e.. Proxy (y) is 
the unique minimizer of the function x i— ^ 2ll^~y|P+/(^)/ the sequence to monitor is 
(Proxy Xjt)jt6N- Now assume that A and B are two nonempty closed convex subsets of X. 
The indicator function la of A takes the value on A, and +oo outside A, and analogously 
for ig. We set f = lb and g = l^- It is then clear that 



(49) Proxy = Pb and Prox 



s 



and that (|48|) turns into 

(50) (Vfc G N) Xfc+i = Txfc, where T = Pa{2Pb - Id) + Id -Pb 

Applying this in X (with A = C and B = D) and recalling (|40]) , we obtain (|47)). Viewed 
directly in X, we obtain the iteration 

(51a) = Tx„ where T := Pc(2Po - Id) + Id - Po = M + PPc - ld)(2Po - Id) ^ 
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and we monitor the sequence 

(51b) (J'DXfc)fceN- 

For convergence results, see ||6l Chapter 26], ||30| , and ||44|. 

Remark 3.14 (D-R vs alternating direction method of multipliers (ADMM)) 

ADMM is a very popular method iT2| that can also be adapted to solve feasibility prob- 
lems for two sets. Suppose we wish to find a point in A n B, where A and B are two 
nonempty closed convex subsets of X. Given uq ^ ^ and fco ^ ^/ ADMM generates three 
sequences (flfc)fc>i/ (&fc)fc6N/ and (wfc)fc6N via 

(52) (Vfc G N) flfc+i := Va{\ - w/c), \+\ := PB(flfc+i + ^k), Wfc+i := % + fl/c+i - 

On the other hand, D-R for this problem, with starting point Xq ^ ^/ produces two 
sequences (xfc)^^^ and (i/fc)fceN via 

(53) (Vfc G N) i/fc := Vb^^, x^+i := PA(2yfc - x^) + Xjt - J/fc- 
Now assume that 

(54) xq = fco £ B and = 0. 

Then i/o = fe^^o = ^0/ = ^^(21/0 - ^o) + - !/o = ^a^o = ^^(^0 - 0) = Pa(&o - "o) = 
a\ = ai + uq and yi = PgXi = Pb{^i + "o) = fci- Furthermore, assume that for some 
k > 1, we have x^ = fl/c + Wfc-i and = ^fc- Then 2i/j;- — x;;- = 2&;c — (wfc-i + flfc) = 
bk — (ik + bk — wjt-i = — + ^/c ~ ^k-i = bk — W/c and x^ — i/jt = ^^fc-i + ^k~bk = u^', 
in turn, this implies Xj^+i = PAi'^Vk ~ ^k) + ^/c ~ Vk = PA{bk ~ + Wfc = '^fc+l + ^^fc and 
i/fc+i = -PfiXit+i = PB(flfc+i + Wfc) = fcfc+i- It follows inductively thal|j 

(55) (^ic)fc>i = (% + "fc-i)/c>i and (i/fc)fc>i = (&/c)fc>i- 
In this sense, 

(56) D-R and ADMM are equivalent. 

See also |[12ll26ll32ll33B for further information on ADMM and related methods. Further- 
more, if A and B are linear subspaces, then (x/Jj^eN/ the sequence governing DR, has the 
update rule 

(57) Xk+^ = iPAPB + PA^PB^)^k 

Since Pa{2Pb - Id) + Id -Pb = Pa{2Pb - Pg - Pg^) + Pr^ = Pa{Pb - Pb^) + PaPb^ + 
P^±Pg± = PaPb + Pa^Pb^- Using, e.g., |[8, Corollary 3.9], one can further show that 
Fix(PAPB + Pa±Pb±) = (A n B) + (a^ n B^). 

^ In fact, this argument works much more generally when the projectors Pa and Pg are replaced by 
arbitrary resolvents and bo is a fixed point of the second resolvent. 
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3.4 Summary of feasibility algorithms 



Name 


Acronym 


Formula 


Monitor 


Cyclic Projections 


CycP 






Cyclic Intrepid Projections 


CycP+ 


(421 




Parallel Projections 


ParP 


m 


ixk)keN 


String-averaging Projections 


SaP 


dSl) 


{Xk)keK 


Extrapolated Parallel Projections 


ExParP 


(SHI) 


{Xk)keN 


Extrapolated Alternating Projections 


ExAltP 


dH]) 




Douglas-Rachford 


D-R 


m 


{^k)keK 



Note that all algorithms in this table proceed by iterating an operator, and monitoring 
either the iterates directly, or some simple version thereof. This is a key point when we 
revisit these algorithms in the next section. 



3.5 Superiorization: between feasibility and best approximation 



The algorithms considered so far are designed to solve the feasibility problem (|34|) . Let 
V G X \ C. We shall discuss algorithms for finding PqV in Section |4] below. This best 
approximation problem is equivalent tc@ solving the optimization problem 

(58) minillx — ylP. 

xec ^ 

The new paradigm of superiorization (see, e.g., fT8\) lies between feasibility and this best 
approximation problem. It is not quite trying to solve (|58|) ; rather, the objective is to find a 
feasible point that is a superior to one returned by a feasibility algorithm. To explain this 
in detail, we assume that T: X ^ X satisfies 

(59) Fix T = C; 
hence, (|58|) is equivalent to 

111 \\ 

(60) min i — ^ • 



Applying the superiorization approach to (|60|), we obtain the following abstract algo- 
rithm: 



We work here with a squared version of ^ because the objective function is then differentiable. 
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(61) 



Algorithm: Superiorization of T 


Data: v eX,£>0 


Result: 




Xq ^ V 




while d(xjt) > e do 




if Xjt — v\\ > then 






X ^ xjt + 9{x]^ — v)/\\x]( — v\\ 




else 






X ^ Xjt 




end 




e ^ 9/2 




if X — v\\ < Xfc — v\\ and d(Tx) < d{x]^) then 






Xfc+i ^ Tx 




end 




k 


^k + 1 


end 





Note that rf: X — > R_|- is a performance function satisfying d{x) = if and only if x G C. 
(In Section[6l we use (|112|) .) With the exception of D-R, each algorithm in Section |3!4l has a 
superiorized counterpart. (It is not clear how D-R should be superiorized because the fixed 
point set of the operator T governing the iteration is generally different from C, the set of 
interest.) We denote the superiorized counterpart of CycP by sCycP, and analogously for 
the other algorithms. For remarks on the numerical performance of these algorithms, see 
Section [631 



4 Best approximation algorithms 
4.1 The problem and notation 

We continue to assume that m G {2, 3, . . .} and that Ci, . . . , C„j are closed convex subsets 
of X such that 

(62) C := Ci n • • • n C,« ^ 0. 
Let u G X. We wish to determine 

(63) Pcv. 
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Before we present a selection of pertinent best approximation algorithms, let us fix some 
notation for ease of use. It is notationally convenient to occasionally work with cyclic 
remainders 

(64) [l]m = l,[2]m = l,...,[m]m = m, [m + l]„, = 1,..., 
taken from {1, . . . , m}. We also require the operator Q defined by |^ 

(65) Q: XxXxX ^ X 

' z, if |0 = and > 0; 

(x,y,z) ^ i ^ + (l + f)(2-y), if^? > Oand^v > p; 



^y + ^(/^(^-y) +f^(z-y))' '^^9 >Oand^v <p, 



where ~ z)/ }i '■= ||x — ylP, v '■= \\y — z|p, and p := }iv — x^- An analogous 

formula holds for Q : X x X x X ^ X. 



4.2 A catalogue of best approximation methods 

As in Section |3l, we present a list of best approximation methods based on projectors and 
comments. Unless stated otherwise, each of these methods produces a main/ governing 
sequence that is converging to PqV. 

Algorithm 4.1 (Halpern-Wittmann (H-W)) Set xq = v and update 

(66) (Vfc G N) x,+i := j^v + t^Pc„.Pc,„^, ■ ■ ■ ^c.^k- 

Remark 4.2 (H-W) This algorithm was introduced by Halpern |l39l while Wittmann ||49l 
proved convergence for the choice of parameters presented in Algorithm 14.11 Many vari- 
ants have been proposed and studied. 

Algorithm 4.3 (Cyclic Dykstra algorithm (CycDyk)) Set xq '•= ^/ cf-^m-i) •= ^-{m-i) •= 

• • • := q-\ := qo := 0, and update 

(67) (Vfc G N) Xk+i := J'C[;,+i]^ {xk + qk+l-m) and q^+l ■= ^fc + qk+l-m - ^fc+l- 

Remark 4.4 (CycDyk) For convergence proofs, see, e.g., |[T4| or HI Theorem 29.2]. See 
also [25} Section 3] for connections to the forward-backward method. 

The following method operates in X. 



If |D = and X < 0, then the output of Q is undefined — this corresponds to the case when C 
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Algorithm 4.5 (Parallel Dykstra algorithm (ParDyk)) Set (yo, zq) = {v, . . . ,v,0, . . . ,0) G 
X X X, and let (y^, zj^) = {yk,\, ■ ■ ■ , yk,m> ^k,\r ■ ■ ■ > ^k,m) ^ X x X be given. Then the next 
iterate is {yk+i,Zk+\) = (j/fc+U'- • •'yfc+i,m/Zfc+U' - • •'Zfc+i,m)/ where the sequence {xk)keN 
to monitor is 



and the update formulas are 

(68b) (Vz G {1, . . . , m}) y^+i^i := Pc,{zk,i + and z^+i,/ := z^^, + % - yk+\,i- 

Remark 4.6 (ParDyk) ParDyk is CycDyk applied to the subsets C and D of X; see, e.g., 
Il2> Theorem 6.1]. 

Remark 4.7 (Dykstra vs ADMM) Let A and B be two nonempty closed convex subsets 
of X, and let y G X. CycDyk for finding P^pB^^^ generates sequences (fljt)fc>i/ (^fc)fc6N/ 
(Pfc))c6N/ and {cik)ke¥i as follow^ Set feg ^/ Po •= •= 0/ and for every G N, update 

(69a) flfc+i := P^l^fc + Vk), Pk+i ■= h + Pk- 

(69b) bk+i := PB(flfc+i + '?/c)/ := %+i + - h+i- 

On the other hand, given G X and yo G X, ADMM for finding a point in A n B — not 
necessarily P^^nB^ — generates three sequences (xfc)/(;>i, {yk)ket<S' and (u/c)a;6N via 

(70) (V/c G N) Xfc+i := PAiyk - Uk), yk+i ■= Psixk+i + u^), Wfc+i := + ^k+i - yk+i- 
Now let us assume that 

(71) B is a linear subspace 

and that uq G B-*-. Then it is well known that the Dykstra update ^69} simplifies to 

(72) flfc+i := PAih + Vk)r h+1 ■= -Pb(«/c+i)/ Pk+i ■= Pk + h- 

because the sequence {cjk)keN lies in B-*- and thus becomes "invisible" when computing 
{bk)keK due to the linearity of Pg. Turning to ADMM, we observe that {u]()keN lies in B-*- 
which simplifies the update for y^+i to y^+i := Psixk+i)- Setting {vk)ken '■= -(i^fc)/c€N/ 
we see that (|70]) turns into 

(73) (Vfc G N) Xfc+i := PAiyk + Vk), yk+i ■= Psixk+i), ^fc+i := Vk + yk+i - x^+i. 

^ Here the sequence {qk)ker^ from (|67t is split into subsequences corresponding to odd and even terms 
for easier readability. 



(68a) 




i=l 
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Comparing (|72|) to ((73|), we see that the update formulas look almost identical: The only 



difference lies in the update formulas of the auxiliary sequences {pk)keN ^rid (^fc)fc6N 
the former works with while the latter incorporates immediately the more recent up- 
date i/fc+i- However, the resulting sequences and hence algorithms appear to be different^ 
Let us now further specialize by additionally assuming that 

(74) A is a linear subspace. 

Then the Dykstra update (|72l) does not require the sequence {pk)keN anymore (because it 
lies in A-*- and it thus plays no role in the generating of (flfc)fc>l)/ and it further simplifies 
to 

(75) flfc+i := PA{bk)> h+i ■= PB{ak+i)- 

Hence, as is well known, Dykstra turns into CycP, which is also known as von Neumann's 
alternating projection method in this special case. On the other hand, these additional as- 
sumptions do not seem to simplify (|73l) . In fact, the behaviour of Dykstra can starkly 
differ from ADMM even in this setting: Suppose that X = R^, that A = R • (1, 1) and 
B = R X {0}. Then the sequence (&A:)fc6N with starting point fco '■= ^ '■= (1/0) turns out 
to be (2~'^,0)fce]N- In contrast (see also Footnote O, we compute the following ADMM 
updates, where i/o '■= ^ := (1/0) and vq := (0,0): Xi = PAiyo + ^o) = PaO->0) = 5(1,1), 

yi = Pb{xi) = i(i,o),ui = uo + yi-xi = |(0,-1),X2 = PA{yi + vi) = P^(i(i,-i)) = 

(0,0), and 1/2 = ^"8(^2) = (0,0). Since X2 = y2 ^ ^ H B, the algorithm terminates when- 
ever feasibility is the implemented stopping criterion. 

Algorithm 4.8 (Haugazeau's algorithm with cyclic projections (hCycP)) Set xq = v, 

and updat^ 

(76) (VfcGN) Xfc+i:=Q(xo,Xfc,Pc[,+^,^Xfc)- 

Remark 4.9 (hCycP) For convergence proofs, see ||40'] or fF, Corollary 29.8]. The general 
pattern for methods that undergo a Haugazeau-type modification is (see [i5J for details) 

(77) (VfcGN) Xfc+i := Q(xo,Xfc,%)Xfc), 

where {Ti)i^j is a family of operators and i: IN ^ J selects which operator is drawn 
at iteration k. This gives rise to many variants; in the following, we shall focus on a 
representative selection. 



It is stated on [12. page 34f] that (in our notation) " ((73|l is exactly Dykstra's alternating projections 
method . . . which is far more efficient than the classical method [of alternating projections] that does not 
use the dual variable v.". This statement appears to be at least ambiguous because the crucial starting 
points are not specified. 
^ Recall llSil- 
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Here is a Haugazeau-type modification of ParP (Algorithm \3A} : 
Algorithm 4.10 (hParP) Set xq := v and update 

(78) (VfcGN) Xk+^:=Qixo,Xk,Txk), where T = 1 (Pc, + ••• + PcJ. 

The following Haugazeau-type modification of D-R operates in the product space X. 

Algorithm 4.11 (hD-R) Let T : X ^ X be the operator governing D-R (see (|51a|) ), and set 
Xq = {v,. . .,v) G X. Given G N and xj^ = (xjt^, . . . , Xj^^m) ^ \ update the governing 
sequence by 

(79a) := Q(xo,Xfc,Txfc), 

and monitor 



^ m 

(79b) xt:= -Y^Xk^i. 



m 



Remark 4.12 (hD-R) To show that x^ Pc{v), use |1 Example 8.2]. 

The next algorithm is a variant of D-R tailored for best approximation. It also operators 
in the product space X. 

Algorithm 4.13 (baD-R) Set xq = {v, . . . ,v) ^ X. Given ?c G N and x^ = (x^^i, . . . , X]^„j) G 
X, update to x^+i = . . . , x^+i^^), where 



^ m 

(80a) -Yl ^^'i 



m 

1=1 



and 

(80b) (Vz G {1, . . . , m}) Xk+ij = Xfc,- - Xfc + Pc.{{v + Ix^ - x^,) /2) . 

4.3 Remarks on baD-R 

We comment on various aspects of baD-R and start with its genesis. 

Remark 4.14 (baD-R) Let A and B be nonempty closed convex subsets of X, and let 
u G X. Revisit Remark |3.13[ in which we showed how the feasibility version of D-R, 
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Algorithm 13.121 arises from the general problem of minimizing f + g,hy choosing f = Lb 
and g = La-^o explain Algorithm I4.13[ we keep f = lb for which 

(81) Proxy = Pg. 

However, this time we take g: x i— > — yjp + fy^(x). Then the (unique) minimizer of 
f + g is indeed Py^nB^- Let x and y be in X. Then y = Prox^x x G (Id +9g)(i/) = 
2y-v + NA{y) ^^{x + v) ey + ^N^(y) =y + N^fy) and so y = PAiji^ + v)); that is, 

(82) Prox^: X ^ P^(i(x + z;)). 

In view of dHD and (|82l) , the general update formula for the Douglas-Rachford splitting 
method, (|48ll , becomes 



(83) (Vfc G N) Xfc+i = Txfc, where T: x^P^((2Pbx-x + u)/2)+x-Pbx. 
Applying this in X (with A = C and B = D) and recalling (|40l) , we obtain ((80l) . 

Remark 4.15 (baD-R = CycDyk = CycP in the doubly affine case) Let us investigate (|83]) 
with starting point u G B further. First, we rewrite it as 

(84) xo:=v e B, and (Vfc G N) y^ := PbXj,, x^+i := x^ - + PA(yfc + (v - Xfc)/2). 
We now assume additionally that 

(85) A and B are linear subspaces of X. 
We claim that for every k > 1, 

(86) Xfc = PaV - PbPav + PaPbPaV T • • • + PAiPsPA^'^ 
which implies 

(87) y,, = {PbPa^v 

because the differences in (|86|) all lie in B-*-. Observe thatyo = PbXq = Pb'v = v since z; G B. 
It follows that xi = xq — yo + ^^(yo + — ^o)/2) = v — v + Pa{v + {v — v) / 0) = Pav 
which shows that (|86l) holds when = 1. Now assume that (|86]) holds for some ?c > 1. 
Then (|87)| also holds. Furthermore, 

(88) Xk-yk = Pav - PbPav ± • • • - {PBPAfv 
and 

(89) z; - Xfc = (Id -P^)z; + (Id -Pa)PbPav + • • • (Id -Pa)(PbPa)''''^^ e A^; 

the latter identity and ((8Zl» yield PA(y/c + (z^ - x^)/l) = PAyk = Pa{PbPaYv. This verifies 
(|86)) with replaced by A: + 1; therefore, by induction, (|86l) and hence (|87l) hold true for 
all A: > 1. In view of (|87l) and (|75l) , we observe that baD-R = CycP = CycDyk in sense 
that the sequences that arise after projection onto B are all identical. Finally, a translation 
argument reduces the doubly affine case to the doubly linear case. 
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Remark 4.16 (baD-R 7^ CycDyk in general) Suppose that A is a nonempty closed con- 
vex subset of X and that B is a subspace of X. Starting both CycDyk and baD-R at v, we 
obtain for G N the update rules 

(90) bo := v,po := 0,%+i := PA{bk + Pk)>Pk+i ■= h + Pk- := -Pb%+i 
and 

(91) xo := v,yk := PsXfc/^fc+i := ^fc - J/fc + PAiVk + {v - Xk)/2), 
respectively. One computes that fli = P^v = x\, p\ = v — PaV, and 

(92) bi = PbPaV = yi. 

It thus follows that ^2 = Pa{PbPa^ + v - Pav) and X2 = Pb^PaV + Pa{PbPaV + (^ - 
Pav)/2)). Consequently, 

(93) b2 = PBPA{PBPAV + v-PAv) and y2 = PBPA{PBPAV+{v-PAv)/2)). 

These two vectors certainly appear to be different — let us now exhibit a simple example 
where they actually are different. We work in the Euclidean plane and thus assume that 
X = R^. Set A := {(^,7) | + (7 — 1)^ ^ l}/ which has the projection formula 

f(^,'7), if^2 + (//-l)2<l; 

^ ^ ^ I (0,1) + =, otherwise, 

and set B := R x {0}. Then 

(95) Pb:R2^r2:(^,,/)^(^,0). 

Set V = (1,0). Then (gH) turns into bi = yi = (\/2/2,0) while (gS)) becomes 

(96) b2 = I , ^ ^ ,0] ^ (0.61180,0) 7^ (0.59718,0) = f 1— ^±i=,0 ) = y2. 

VV22-8\/2 ) \}^\\-l^fl ) 

Let us mention in passing that similar examples exist in the product space setting of Sec- 
tion |3]1] when some of the constraint sets are balls. 
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4,4 Summary of best approximation algorithms 



Name 


Acronym 


Formula 


Monitor 


Halpem-Wittmann 


H-W 


(|66D 




Cyclic Dykstra 


CycDyk 


(67) 




Parallel Dykstra 


ParDyk 


(68) 


iXk)keN 


Haugazeau-like Cyclic Projections 


hCycP 


(76) 


{Xk)keN 


Haugazeau-like Parallel Projections 


hParP 


dZSl) 


{Xk)keN 


Haugazeau-like Douglas-Rachford 


hD-R 


(|79D 


{Xk)keN 


best approximation with Douglas-Rachford 


baD-R 


(80) 


{^k)keN 



5 Nonconvexity 

In this section, we provide some remarks on the possible absence of convexity. 



5.1 Nonconvex slope constraints 

We now consider a case with nonconvex constraints. This type of constraints does occur 
frequently in applications; however, the body of convergence results is sparse and, to the 
best of our knowledge, all results are local, i.e., convergence of the iterates is guaranteed 
only when the starting point is already sufficiently close to the set of solutions. (See, e.g., 
IITOl and the references therein.) We use the convex constraints from Section |2] with one 
crucial modification: The slope constraints (|T3l) is tightened by additionally imposing a 
minimum absolute value slope, i.e., we assume the existence of two vectors a = (a;) and 
b = (|6f) in R^r^ such that 

(97) (Vz G I) Si := {x = (xi,. . .,x„) G R" | at > \x^+l - x,\ > ^6,}. 



Note that if j6j > 0, then Sj is not convex. Analogously to (|T8|) , we aggregate these sets to 
obtain the general nonconvex constraint set 

(98) S= fl Si. 

ie{l,...,n-l} 
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Let i G {1, ... ,n — 1} and x G X. Arguing as in Section |2.3.1[ we obtain the counterpart 
of (|T6l) and see that Ps^x is different from x at most in positions i and i + 1, anc^ 



X/+1 - Xi < di 




Figure 2: General slope constraint in the nonconvex case. 



if ^ Xi ^i'f 

if Xi - di < Xi+i < Xi - ^i- 
if Xi - ^i < Xi+i < Xi) 
if Xi < Xij^i < Xi + jS;-; 
if Xi + ^i < Xi^i < Xi + di) 
if Xi + di < Xf+i. 

This formula allows us to deal with the general case as in Section 12.3. 2t thus, we obtain 
the counterpart of Proposition 12. 2[ 

^ The astute reader will note that when /5/ > and Xj = then 

is not single-valued; indeed, in view of the nonconvexity of S, and the famous Bunt-Motzkin Theorem (see, 
e.g., [6. Corollary 21.13]), this is to be expected. For the actual implementation, one chooses an arbitrary 
element from this set. 



(99) 



((Ps,x).(Ps,x)m) 



\{xi + x,+i + ^i, Xi + x;+i - j6j), 

^(X; + X;_|_X — /3j, X; + Xf + 1 + /3;), 

(Xf, Xj^x)/ 
, ^(Xj + X;_|_x — di, Xi + Xj+1 + di), 
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Let us also record the intrepid counterpart of (|99|): 
(100) 

l{xi + Xi+i + lioci + ^i),Xi + - \{DCi + j6,)), if < Xi + (|6j - 3a,)/2; 

+ a;;,Xj - cii), if + (|6/ - 3a:/) /2 < < - a,; 
(x/,x,+i), if X/ - DCi < X/+1 < X, - j6/; 

(x/+i + j6/, Xi - ^i), if X/ - < Xi+i < Xi + miri{0, (a, - 3^i)/2}; 

l{xi + Xi+i + \{oci + p>i),Xi + x,+i - \{DCi + |6,)), if Xi + (a; - 3|6,)/2 < x,+i < x;; 

< 

\{xi + x,-+i - i(a/ + ^i),Xi + x,+i + \{DCi + |6/)), if X; < X/+1 < Xi + (3|6/ - a;)/2; 
(xj+i - ^i,Xi + |6,), if X; + max{0, (3|6, - a;/)/2} < x,+i < x, + ^i; 
(x/,x/+i), if x, + ^i < Xi+i < Xi + (Xi) 
{xi+i - Oii,Xi + oii), if Xi + oii < X/+1 < Xi + (3a;, - j6/)/2; 
J(x, + Xi+i - j{(Xi + p>i),Xi + x,+i + i(a, + p>i)), if X, + (3a/ - j6,)/2 < X/+i. 

There are at least 8 cases; however, the two cases in the middle arise if and only if 3/3, > a;,. 

Having recorded all required formulas, one can now experiment with the performance 
of the corresponding algorithms. In the physics community, Veit Elser has championed 
especially D-R with great success (see, e.g., | 31il and Il38ll ). Because of its central impor- 
tance, we consider in the following subsection a curious cycling behaviour of D-R. 

5.2 Representation of nonconvex constraints and cycling for D-R 

Suppose momentarily that Ci, C2, C3 are nonempty closed convex subsets of X such that 
C = Ci n C2 n C3 7^ 0. Assume further that Ci n C2 is still simple enough to admit a 
formula for the projector PcinC2- The projection methods considered find points in C; it 
does not matter if we work with either the original three sets C\, C2, C3 or with C\ n C2 
and C3. This situation changes dramatically if we are dealing with nonconvex constraints. 
Performance of projection methods crucially depends on the representation of the con- 
straints. 

In fact, the example developed next was discovered numerically. It forced us us to deal 
with nonconvex constraints differently and led to the formula (|99l) . 
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Example 5.1 (cycling for D-R) Suppose that X = R^, that a > ^ > 0, and set 

(101) Ci := {{xi,X2) I \xi — X2I < i^} and C2 := {{xi,X2) \ \xi — X2I > j^}- 
Clearly, Ci n C2 7^ 0; in fact, (|99)) allows us to record a formula for PcinC2- Set 

(102) e := ^6/5, 
take an arbitrary ^ G IR, and set 

(103) xo = (X04, xo,2) = (^, ^ - £, ^ - 2e, ^ + e) G r2 X R^. 

Then the sequences (xfc)fceN and (xjt)fce]N generated by (|47|) do not converge; indeed, they 
cycle as follows: 

(104) (VfcGN) X2fc = (^,^-e,^-2e,^ + £), X2fc+i = - e,^,^ + e,^ - 2e) 
and 

(105) (VfcGN) X2fc = (^-e,a %m = (^^^-e)- 

Proof. It is clear that xq = — e, ^). Observe that 

U^^nl if|^-7| >/3; 

(106) (V(^,//)gR2) Pc^(^,,y) = J i(^ + ;y-/3,^ + ,/ + /3), if l^-i/l </3and^<//; 

[l(^ + // + /3,^ + //-/3), if |^-//| </3and^>//. 

By definition, 

(107a) X14 = xo,i - ^0 + Pci (2x0 - ^0,1) 

(107b) = (^, ^ - e) - (^ - £, a + Pc, (2(^ - £, a - (^, ^ - e) ) 

(107c) = (£,-£)+PQ(^-2e,^ + e). 

Now I (^ - 2e) - (^ + e) I = 3e = (3/5)j6 < |6 < a, so Pq does not modify (^ - 2e, ^ + e) 
and we conclude that 

(108) xi,i = (£, -£) + (^ - 2e, ^ + £) = (^ - £, ^). 

Next, 

(109a) xi,2 = xo,2 - + J'cj (2xo - ^0,2) 

(109b) = (^ - 2e, ^ + e) - (^ - e, ^) + Pc, (2(^ - e, ^) - - 2e, ^ + e) ) 

(109c) =(-£,£) +Pc,(^,^-e). 
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Now 1^ - - e) I = e = j6/5 < |6, so (|T06ll yields 



(110) Pc^i^,^- e) = \{n - e + /3,2^ - £ - /3) = (^ + 2£,^ - Se). 
It follows that xi,2 = (-£, e) + + 2e, ^ - 3e) = + e, ^ - 2e). Altogether, 

(111) xi = (xi4,xi,2) = (^-e,^,^ + e,^-2e) and Xi = (^,^-e). 
Arguing similarly, we obtain that X2 = xq. 



Remark 5.2 If D-R is started at different points, as we recommend in Algorithm |3.12[ 
then the iterates eventually settle into this cycling behaviour described in Example l5.ll 



6 Numerical results 

6.1 Experimental setup and stopping criteria 

6.1.1 Generating the set of test problems 

We randomly generate 100 test problems, namely road splines in groups centered 
around length L G {0.5,1,5,10,20} (unit km). For each group, we select a de- 
sign speed V G {30,50,80,100} (with unit km/h) and maximum elevations ^max ^ 
{30, 60, 100, 120, 150} (unit m). We then generate spline points { {t\, v\),. . . , {tn, Un) } such 
thatu G [0,^max]" and (Vz G {l,...,n-l}) \\{t„v,) - {t,^x,Vr^\)\\ > 0. 625 V, where 
n G [L/(3min{0.625V,30}),...,l + L/(1.5min{0.625V,30})] nN. The resulting splines 
correspond to rather challenging road profiles and are therefore ideal for testing. 

6.1.2 Stopping criteria 

The constraint sets Cf are as in Section |2] (and Section 15.11 in the nonconvex case). Let 
{xk)keiN be a sequence to monitor generated by an algorithm under consideration and set 
e := 5 • lO""^. The feasibility performance measure we employ is 



Note that d returns the value precisely at points drawn from C = Ci n • • • n C^. If we are 
interested in feasibility only, then we terminate when 



(112) 




(113) 



d{xk) < e; 
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in case of best approximation {Pc'u, where y G X is given), we require additionally that 
— II < e. We cap the maximum number of iterations at kj^^^ — 5000. 



6.2 Generation of plots and tables 

Let V be the set of test problems and A be the set of algorithms. Let (x['''''^)jte]N be the 
monitored sequence generated by algorithm a ^ A applied to the problem p ^ V. 



6.2.1 Performance plots 



To compare the performance of the algorithms, we use performance profiles: for every a ^ A 
and for every p ^ V, we set 



(114) ra,p := . T ,, ^ > 1' 



min {kiji p \ a' ^ A} 



where ka,p G {1,2, . . .,fcmax} is the number of iterations that a requires to solve p (see 
Section [6.1. 2[) . If ra,p = 1, then a uses the least number of iterations to solve problem p. If 
i'a,p > 1/ then a requires ra,p times more iterations for p than the algorithm that uses the 
least iterations for p. For each algorithm a ^ A, we plot the function 

card \p eV \ log^fr^ „) < k) 
(115) p. : IR+ ^ [0, 1] : K ^ card^ ' 

where "card" denotes the cardinality of a set. Thus, pa (k) is the percentage of problems 
that algorithm a solves within factor 2^ of the best algorithms. Therefore, an algorithm 
a ^ A is "fast" if pai^) is large for k small; and a is "robust" if pai^) is large for k large. 
For further information on performance profiles we refer the reader to [1291 . 



6.2.2 Runtime plots 

The sequence (rf(x['''''''))fc6N measures the runtime progress of a on p with respect to the 
feasibility performance measure d (see (|112|) ). To get a sense of the average (logarithma- 
tized) progress of each algorithm a ^ Aat iteration k, we follow [[23l and plot the values 
of relative proximity function, which is defined by 

(116a) N ^ R: ^ lOlog^o ^ ^'(4"''^)) 

peV 
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(116b) 



10 log 



10 



—y 

card V ^ 



{a,p) {a,p) 











6.2.3 Distance tables 

For each algorithm a ^ A and problem p ^ V, assume that termination occurred at 
iteration k(a, p). We compute the normalized distance to v by 



(117) 



\v — X 



{a,v) I 



max 



X 



i' G A], 



otherwise. 



This allows us to consider the collection of normalized distances {^(^a,p))peV for each al- 
gorithm a ^ A. Statistical values are recorded for each algorithm in a table allowing us to 
compare best approximation performance. 



6.3 Results for feasibility algorithms 



6.3.1 The convex case 



In this section, we record the results for the convex setting of Section |2] using the algo- 
rithms of Section 1231 
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Figure 3: Performance profiles for feasibility algorithms in the convex case (see Sec- 
tion [6231) 
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Figure 4: Runtime plots for feasibility algorithms in the convex case (see Section [6.2.2[) 



Algo. 


A Min. 


Al^' Qrt. 


A Median 


A3'-^ Qrt. 


A Max. 


A Mean 


A Std.dev. 


CycP 


0.0563 


0.3039 


0.3696 


0.4806 


0.9301 


0.3993 


0.1568 


CycP+ 


0.0563 


0.2862 


0.3589 


0.4403 


0.7707 


0.3758 


0.1399 


ParP 


0.0517 


0.2973 


0.3687 


0.4770 


0.9191 


0.3990 


0.1612 


SaP 


0.0563 


0.3033 


0.3693 


0.4795 


0.9292 


0.3993 


0.1577 


ExParP 


0.0522 


0.3012 


0.3697 


0.4767 


0.9257 


0.3981 


0.1572 


ExAltP 


0.0522 


0.3012 


0.3697 


0.4767 


0.9257 


0.3973 


0.1558 


D-R 


0.0721 


0.3833 


0.4555 


0.5486 


1.1353 


0.4807 


0.1861 



Table 1: Statistical data for feasibility algorithms in the convex case (see Section [6. 2.3|) 



6.3.2 The nonconvex case 

The plots and tables are analogous to those of Section [6 .3 . 1 1 except that we work with the 
nonconvex slope constraint (see Section |5J|). 
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Figure 5: Performance profiles for feasibility algorithms in the nonconvex case (see Sec- 
tion [6231) 
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Figure 6: Runtime plots for feasibility algorithms in the nonconvex case (see Section |6.2.2|) 
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A Min 




A MpHipin 




A Mpiv 


A Mppin 


A StH Hpv 




U.UJUD 




n '^871 

U.OO/ i 








U.iuZo 


CycP+ 


0.0566 


0.2520 


0.3182 


0.3708 


0.5873 


0.3135 


0.1047 


ParP 


0.0522 


0.3004 


0.3919 


0.4853 


0.9275 


0.4137 


0.1653 


SaP 


0.0566 


0.3036 


0.3820 


0.4824 


0.9292 


0.4027 


0.1574 


ExParP 


0.0527 


0.3031 


0.3796 


0.4786 


0.9306 


0.4026 


0.1603 


ExAltP 


0.0526 


0.3022 


0.3698 


0.4771 


0.9243 


0.3981 


0.1566 


D-R 


0.0756 


0.3835 


0.4567 


0.5486 


1.1044 


0.4810 


0.1846 



Table 2: Statistical data for feasibility algorithms in the nonconvex case (see Section [6.2.3)1 



6.3.3 Conclusions 

Overall, CycP+, ExAltP, and D-R emerge as good algorithms for feasibility. CycP+ yields 
solutions closest to the initial spline and is among the most robust. In the convex case, 
CycP+ is also the fastest algorithm; in the nonconvex case, D-R is faster. Both ExAltP and 
D-R are good choices when the number of iterations is small. 

6.4 Results for superiorized feasibility algorithms 
6.4.1 The convex case 

In this section, we record the results for the convex setting of Section |2] using superiorized 
feasibility algorithms (see Section 133]) . 
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Figure 7: Performance profiles for superiorized feasibility algorithms in the convex case 
(see Section [6.2. 1|) 
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Algo. 


A Min. 


Al^' Qrt. 


A Median 


A3"^ Qrt. 


A Max. 


A Mean 


A Std.dev. 


sCycP 


0.0519 


0.3039 


0.3696 


0.4805 


0.9301 


0.3992 


0.1570 


sCycP+ 


0.0563 


0.2863 


0.3599 


0.4410 


0.7707 


0.3761 


0.1399 


sParP 


0.0517 


0.2973 


0.3659 


0.4762 


0.9191 


0.3940 


0.1559 


sExParP 


0.0563 


0.3039 


0.3711 


0.4825 


0.9301 


0.4018 


0.1559 


sExAltP 


0.0519 


0.3039 


0.3711 


0.4825 


0.9301 


0.4017 


0.1561 


sSaP 


0.0519 


0.3033 


0.3693 


0.4795 


0.9292 


0.3985 


0.1567 



Table 3: Statistical data for superiorized feasibility algorithms in the convex case (see 
Section [6^ 



6.4.2 The nonconvex case 



The plots and tables are analogous to those of Section [6.4. 1 1 except that we work with the 
nonconvex slope constraint (see Section |5l^. 




Figure 9: Performance profiles for superiorized feasibility algorithms in the nonconvex 
case (see Section [6.2. 1|) 
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Figure 10: Runtime plots for superiorized feasibility algorithms in the nonconvex case 
(see Section [6.2.2[) 



Algo. 


A Min. 


Al^' Qrt. 


A Median 


A3"^ Qrt. 


A Max. 


A Mean 


A Std.dev. 


sCycP 


0.0566 


0.3039 


0.3754 


0.4759 


0.9301 


0.4002 


0.1563 


sCycP+ 


0.0566 


0.2885 


0.3594 


0.4388 


0.7708 


0.3765 


0.1392 


sParP 


0.0522 


0.3004 


0.3754 


0.4744 


0.9275 


0.3973 


0.1554 


sExAltP 


0.0523 


0.3042 


0.3754 


0.4762 


0.9301 


0.4030 


0.1555 


sExParP 


0.0566 


0.3042 


0.3754 


0.4762 


0.9301 


0.4030 


0.1554 


sSap 


0.0522 


0.3036 


0.3714 


0.4734 


0.9292 


0.3978 


0.1556 



Table 4: Statistical data for superiorized feasibility algorithms in the nonconvex case (see 
Section [6:Z3l) 



6.4.3 Conclusions 

We clearly see that sCycP+ is not only the fastest but also the most robust superiorized 
feasibility algorithm. 
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6.5 Results for best approximation algorithms 
6.5.1 The convex case 

In this section, we record the results for the convex setting of Section |2] using the best 
approximation algorithms of Section HI 
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Figure 11: Performance profiles for best approximation algorithms in the convex case (see 
Section [623]) 
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A] cm 


A Min 




A MpHipin 




A Mpiv 


A Mppin 


A StH Hpv 


in— vv 


U.UDAD 


n 9Q78 




U.^OUU 




U.oyDy 


n 1 R7n 

U.LD/V 


CycDyk 


0.0519 


0.3022 


0.3695 


0.4862 


0.9379 


0.4006 


0.1584 


ParDyk 


0.0517 


0.2962 


0.3656 


0.4799 


0.9379 


0.3972 


0.1576 


hCycP 


0.0519 


0.2986 


0.3645 


0.4801 


0.9251 


0.3962 


0.1569 


hParP 


0.0517 


0.3042 


0.3720 


0.4864 


0.9379 


0.4016 


0.1585 


hD-R 


0.0519 


0.3042 


0.3720 


0.4864 


0.9337 


0.4013 


0.1579 


D-Rba 


0.0517 


0.2962 


0.3656 


0.4799 


0.9379 


0.3972 


0.1576 



Table 5: Statistical data for best approximation algorithms in the convex case (see Sec- 
tion [l^iSl) 



6.5.2 The nonconvex case 



The plots and tables are analogous to those of Section [6 .5 . 1 1 except that we work with the 
nonconvex slope constraint (see Section |5J|) . 
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Figure 13: Performance profiles for best approximation algorithms in the nonconvex case 
(see Section [6.2.1|| 
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Algo. 


A Min. 


Al^' Qrt. 


A Median 


A3'-^ Qrt. 


A Max. 


A Mean 


A Std.dev. 


H-W 


0.0525 


0.3019 


0.3794 


0.4825 


0.9252 


0.3998 


0.1570 


CycDyk 


0.0667 


0.3053 


0.3861 


0.4863 


0.9385 


0.4079 


0.1539 


ParDyk 


0.0667 


0.3006 


0.3846 


0.4863 


0.9385 


0.4060 


0.1529 


hCycP 


0.0522 


0.3022 


0.3782 


0.4808 


0.9253 


0.3998 


0.1568 


hParP 


0.0520 


0.3057 


0.3863 


0.4864 


0.9385 


0.4069 


0.1570 


hD-R 


0.0652 


0.3057 


0.3825 


0.4864 


0.9342 


0.4074 


0.1540 


baD-R 


0.0667 


0.3006 


0.3846 


0.4863 


0.9385 


0.4060 


0.1529 



Table 6: Statistical data for best approximation algorithms in the nonconvex case (see 
Section i^S]) 



6.5.3 ParDyk vs baD-R 

The plots and tables of Sections 16.5.11 and 16.5.21 suggest that ParDyk and baD-R are the 
same algorithms. This is not true for these algorithms in their full generality as we pointed 
out in Remark |4.16[ A possible explanation for this identical performance could be that 
the algorithms produce sequences that stay either outside or inside the halfspaces com- 
prising the constraints. If that is indeed the case, then the projection onto the halfspace 
would be indistinguishable from the projection onto either a hyperplane or the entire 
space; consequently, the iterates would be identical by Remark [4. 151 
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6.5.4 Conclusions 



hCycP is a good robust choice for both convex and nonconvex problems; CycDyk does 
well for convex problems, and H-W does well for nonconvex problems. 

Examining the data for all three types of algorithms (feasibility, superiorization, best 
approximation), we can make a compelling case for CycP+, which emerges as the best 
overall algorithm with respect to speed, robustness, and best approximation properties. 

6.6 Projection methods vs Linear Programming (LP) algorithms 

In Section |3^ we made the case for projection methods and commented on their compet- 
itiveness with other optimization methods. To illustrate this claim, we interpreted ((33)) 
as the constraints of a Linear Programming (LP) problem. We then solved our test prob- 
lems with the GNU Linear Programming Kit (GLPK) [34J. Two objective functions were 
employed: x i— ^ 0, which corresponds to the feasibility problem ((34|) , and x i— ^ ||x — z;||i, 
which is similar to the best approximation problem^ (|35|) . We call these methods GLPKO 
and GLPKl, respectively 1481 page 257]. We shall compare these two methods against the 
overall best projection method, CycP+. 



Name 


Acronym 


Cyclic Intrepid Projections 


CycP+ 


0-objective function LP 


GLPKO 


£x -minimization LP 


GLPKl 



Here, ||x||i = J2^=i denotes the fj-norm of x e X. 
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Figure 15: Performance profiles for convex feasibility problems comparing the number of 
iterations of CycP+, GLPKO and GLPKl until convergence (see Section [6 .2 .11) 



Algo. AMin. Al^' Qrt. A Median A3"' Qrt. A Max. A Mean A Std.de v. 

CycP+ 0.0563 0.2862 0.3589 0.4403 0.7707 0.3758 0.1399 

GLPKO 0.1938 0.6636 0.8030 0.9281 1.0039 0.7767 0.1813 

GLPKl 0.0000 0.3062 0.3813 0.5001 0.9406 0.4076 0.1625 



The nonconvex slope constraints from Section |5H can be handled by using binary vari- 
ables; this leads to a Mixed Integer Linear Programming (MIP) problem (see, e.g., Il43i Chap- 
ter 12]), which GLPK O can handle as well. 



Table 7: Statistical data in the convex case (see Section [6.2.3|) 
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Figure 16: Performance profiles for nonconvex feasibility problems comparing the num- 
ber of iterations of CycP+, GLPKO and GLPKl until convergence (see Section [6.2.1[) 



Algo. 


A Min. 


Al^' Qrt. 


A Median 


A3'-^ Qrt. 


A Max. 


A Mean 


A Std.dev. 


CycP+ 


0.0566 


0.2520 


0.3182 


0.3708 


0.5873 


0.3135 


0.1047 


GLPKO 


0.1938 


0.5371 


0.6245 


0.8240 


6.2070 


0.8226 


0.7464 


GLPKl 


0.0840 


0.3127 


0.3847 


0.5001 


2.1468 


0.4261 


0.2325 



Table 8: Statistical data in the nonconvex case (see Section [6. 2.3[) 



7 Concluding remarks 

Using the practical example of road design, we formulated (convex and nonconvex) fea- 
sibility and best approximation problems. We studied projection methods and imple- 
mented them to solve the feasibility problems. A clear winner emerged, even when com- 
pared to LP solvers: CycP+, the method of cyclic intrepid projections. 

In the future, we plan to study the influence of parameters on the performance of the 
algorithms presented. The design and testing of hybrid methods that aim to combine 
advantageous traits of various algorithms is also of considerable interest. Finally, rigorous 
convergence statements about CycP+ in the nonconvex setting await to be discovered. 
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